library(AnNet)
library(synaptome.db)
library(pander)
library(ggplot2)
library(ggrepel)
scale <- function(x, VALUE=NULL){
x = as.numeric(as.vector(x))
xmin <- min(x,na.rm=T)
xmax <- max(x,na.rm=T)
if( is.null(VALUE) ){
x <- x-xmin
x <- ifelse(!is.na(x), x/(xmax-xmin), NA)
return(x)
}
value = as.numeric(as.vector(value)[1])
value = value-xmin
value = ifelse(!is.na(value), value/(xmax-xmin), NA)
return(value)
}The package allows building the network from the data frame, where the rows correspond to the edges of the graph, or for the list of nodes (genes), for which the information will be retrieved from the SynaptomeDB.
cid<-match('Presynaptic',getCompartments()$Name) # Let's get the ID for presynaptic compartment
cid
#> [1] 2
t<-getAllGenes4Compartment(cid) # Now we need to collect all the gene IDs for presinaptic compartment
dim(t)
#> [1] 2304 8
head(t)
#> # A tibble: 6 × 8
#> GeneID MGI HumanEntrez MouseEntrez RatEntrez HumanName MouseName RatName
#> <int> <chr> <int> <int> <int> <chr> <chr> <chr>
#> 1 1 MGI:1277… 1742 13385 29495 DLG4 Dlg4 Dlg4
#> 2 2 MGI:88256 815 12322 25400 CAMK2A Camk2a Camk2a
#> 3 3 MGI:96568 9118 226180 24503 INA Ina Ina
#> 4 4 MGI:98388 6711 20742 305614 SPTBN1 Sptbn1 Sptbn1
#> 5 5 MGI:88257 816 12323 24245 CAMK2B Camk2b Camk2b
#> 6 6 MGI:1344… 1740 23859 64053 DLG2 Dlg2 Dlg2
gg<-buildFromSynaptomeByEntrez(t$HumanEntrez) # finally, build the graph using respecctive gene EntrezIDs as node IDs
summary(gg)
#> IGRAPH be669ac UN-- 1780 6620 --
#> + attr: name (v/c)
gg<-annotateGeneNames(gg)
#> 'select()' returned 1:1 mapping between keys and columns
summary(gg)
#> IGRAPH be669ac UN-- 1780 6620 --
#> + attr: name (v/c), GeneName (v/c)
head(V(gg))
#> + 6/1780 vertices, named, from be669ac:
#> [1] 273 6455 1759 1785 26052 2923
head(V(gg)$GeneName)
#> [1] "AMPH" "SH3GL1" "DNM1" "DNM2" "DNM3" "PDIA3"alg<- 'louvain'
gg<-calcClustering(gg,alg = alg){r cluster.graph, cache=TRUE}
gg<-calcAllClustering(gg)
gg <- calcCentrality(gg)Build consensus matrix for louvain clustering
conmat<-makeConsensusMatrix(gg,N=100,mask = 10,alg = alg,type = 2)steps <- 100
Fn <- ecdf(conmat[lower.tri(conmat)])
X<-seq(0,1,length.out=steps+1)
cdf<-Fn(X)
dt<-data.frame(cons=X,cdf=cdf)
ggplot(dt,aes(x=cons,y=cdf))+geom_line()+
theme(
axis.title.x=element_text(face="bold",size=rel(2.5)),
axis.title.y=element_text(face="bold",size=rel(2.5)),
legend.title=element_text(face="bold",size=rel(1.5)),
legend.text=element_text(face="bold",size=rel(1.5)),
legend.key=element_blank())+
theme(panel.grid.major = element_line(colour="grey40",size=0.2),
panel.grid.minor = element_line(colour="grey40",size=0.1),
panel.background = element_rect(fill="white"),
panel.border = element_rect(linetype="solid",fill=NA))clrob<-getRobustness(gg,alg,conmat)
pander(clrob)| alg | C | Cn | Crob | CrobScaled |
|---|---|---|---|---|
| louvain | 1 | 189 | 0.09604 | 1 |
| louvain | 2 | 217 | 0.09 | 0.4831 |
| louvain | 3 | 282 | 0.08741 | 0.2615 |
| louvain | 4 | 139 | 0.08849 | 0.3533 |
| louvain | 5 | 121 | 0.09316 | 0.7536 |
| louvain | 6 | 163 | 0.08759 | 0.2769 |
| louvain | 7 | 30 | 0.08436 | 0 |
| louvain | 8 | 225 | 0.08831 | 0.338 |
| louvain | 9 | 150 | 0.0866 | 0.1922 |
| louvain | 10 | 108 | 0.09031 | 0.5097 |
| louvain | 11 | 106 | 0.08519 | 0.0714 |
| louvain | 12 | 19 | 0.09012 | 0.4934 |
| louvain | 13 | 31 | 0.09251 | 0.6978 |
bm<-getBridgeness(gg,alg,conmat)
#> calculating Bridgeness for: louvain
pander(head(bm))| ENTREZ.ID | GENE.NAME | BRIDGENESS.louvain |
|---|---|---|
| 273 | AMPH | 0.0957579145556555 |
| 6455 | SH3GL1 | 0.452106501493934 |
| 1759 | DNM1 | 0.0892319279050531 |
| 1785 | DNM2 | 0.254828113224788 |
| 26052 | DNM3 | 0.147644688926001 |
| 2923 | PDIA3 | 0.548658380347371 |
VIPs=c('8495','22999','8927','8573','26059','8497','27445','8499')
# VIPs=c('81876','10890','51552','5874','5862','11021','54734','5865','5864',
# '9522','192683','10067','10396','9296','527','9114','537','535',
# '528','51382','534','51606','523','80331','114569','127262','57084',
# '57030','388662','6853','6854','8224','9900','9899','9145','9143',
# '6855','132204','6857','127833','6861','529','526','140679','7781',
# '81615','6844','6843')
indx <- match(V(gg)$name,VIPs)
group <- ifelse( is.na(indx), 0,1)
MainDivSize <- 0.8
xmin <- 0
xmax <- 1
ymin <- 0
ymax <- 1
Xlab <- "Semilocal Centrality (SL)"
Ylab <- "Bridgeness (B)"
X <- as.numeric(igraph::get.vertex.attribute(gg,"SL",V(gg)))
X <- scale(X)
Y <- as.numeric(as.vector(bm[,dim(bm)[2]]))
lbls <- ifelse(!is.na(indx),V(gg)$GeneName,"")
dt<-data.frame(X=X,Y=Y,vips=group,entres=V(gg)$name,name=V(gg)$GeneName)
dt_vips<-dt[dt$vips==1,]
dt_res<-dt[dt$vips==0,]
##--- baseColor of genes
baseColor="royalblue2"
##--- SPcolor, colour highlighting any 'specical' genes
SPColor="royalblue2"
##--- PSDColor, colour of core PSD95 interactor genes
PSDColor="magenta"
ggplot(dt,aes(x=X,y=Y,label=name))+#geom_point()+
geom_point(data=dt_vips,
aes(x=X,y=Y),colour=baseColor,size = 7,shape=15,show.legend=F)+
geom_point(data=dt_res,
aes(x=X,y=Y, alpha=(X*Y)), size = 3,shape=16,show.legend=F)+
geom_label_repel(aes(label=as.vector(lbls)),
fontface='bold',color='black',fill='white',box.padding=0.1,
point.padding=NA,label.padding=0.15,segment.color='black',
force=1,size=rel(3.8),show.legend=F,max.overlaps=200)+
labs(x=Xlab,y=Ylab,title=sprintf("%s",alg))+
scale_x_continuous(expand = c(0, 0), limits = c(xmin, xmax)) +
scale_y_continuous(expand = c(0, 0), limits = c(ymin, ymax))+
theme(
axis.title.x=element_text(face="bold",size=rel(2.5)),
axis.title.y=element_text(face="bold",size=rel(2.5)),
legend.title=element_text(face="bold",size=rel(1.5)),
legend.text=element_text(face="bold",size=rel(1.5)),
legend.key=element_blank())+
theme(panel.grid.major = element_line(colour="grey40",size=0.2),
panel.grid.minor = element_line(colour="grey40",size=0.1),
panel.background = element_rect(fill="white"),
panel.border = element_rect(linetype="solid",fill=NA))+
geom_vline(xintercept=0.5,colour="grey40",size=MainDivSize,linetype=2,show.legend=F)+
geom_hline(yintercept=0.5,colour="grey40",size=MainDivSize,linetype=2,show.legend=F)
#> Warning: ggrepel: 1 unlabeled data points (too many overlaps). Consider
#> increasing max.overlapslibrary(plotly)
p<-ggplot(dt,aes(x=X,y=Y,label=name))+#geom_point()+
geom_point(data=dt_vips,
aes(x=X,y=Y),colour=baseColor,shape=15,show.legend=F)+
geom_point(data=dt_res,
aes(x=X,y=Y, alpha=(X*Y)),shape=16,show.legend=F)+
labs(x=Xlab,y=Ylab,title=sprintf("%s",alg))+
scale_x_continuous(expand = c(0, 0), limits = c(xmin, xmax)) +
scale_y_continuous(expand = c(0, 0), limits = c(ymin, ymax))+
theme(
axis.title.x=element_text(face="bold",size=rel(2.5)),
axis.title.y=element_text(face="bold",size=rel(2.5)),
legend.title=element_text(face="bold",size=rel(1.5)),
legend.text=element_text(face="bold",size=rel(1.5)),
legend.key=element_blank())+
theme(panel.grid.major = element_line(colour="grey40",size=0.2),
panel.grid.minor = element_line(colour="grey40",size=0.1),
panel.background = element_rect(fill="white"),
panel.border = element_rect(linetype="solid",fill=NA))+
geom_vline(xintercept=0.5,colour="grey40",size=MainDivSize,linetype=2,show.legend=F)+
geom_hline(yintercept=0.5,colour="grey40",size=MainDivSize,linetype=2,show.legend=F)
ggplotly(p)
#> Warning: `gather_()` was deprecated in tidyr 1.2.0.
#> Please use `gather()` instead.
#> This warning is displayed once every 8 hours.
#> Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.